% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Load Data
FileName = ['AllTradingDate19962021_Monthly.mat'];
load(FileName, ...
     'Target_AllDate');             
clear FileName

for d = 1:length(Target_AllDate)
    Target_Date = Target_AllDate(d, 1);
    Target_TTM = Target_AllDate(d, 3);

    % Load Data
    FileName = ['OP' num2str(fix(Target_Date / 10000)) '_' num2str(fix(rem(Target_Date, 10000) / 100)) '.txt'];
    Data = load(FileName);     
    clear FileName

    Index_ID = 1;
    Index_Date = 2;
    Index_TTM = 3;
    Index_CPFlag = 4;
    Index_K = 5;
    Index_S = 6;
    Index_OP_Bid = 7;
    Index_OP_Ask = 8;
    Index_OI = 9;
    Index_V = 10;
    Index_RF = Index_V + 1;                                                
    Index_DY = Index_V + 2;                                                
    Index_IS = Index_V + 3;                                                
    Index_IS_ADJ = Index_V + 4;                                            
    Index_S_ADJ = Index_V + 5;                                             
    Index_IV_Bid = Index_V + 6;                                            
    Index_IV_Ask = Index_V + 7;                                            

    Date_EXP_WeekDay = weekday(datenum(num2str(Data(:, Index_Date)), 'yyyymmdd') ...
                               + Data(:, Index_TTM));      
    Data(Date_EXP_WeekDay==7, Index_TTM) = Data(Date_EXP_WeekDay==7, Index_TTM) - 1;       
    clear Date_EXP_WeekDay            
    
    Index = find((Data(:, Index_Date)==Target_Date) & ...
                 (Data(:, Index_TTM)==Target_TTM));
    Data = Data(Index, :);
    Data(:, Index_TTM) = Data(:, Index_TTM) - 1;     
    clear Index
 
    % Risk-Free Rate  
    FileName = ['RiskFreeRate19962021.txt'];
    Data_RF = load(FileName);
    clear FileName

    Data(:, Index_RF) = RF_TTM(Data_RF, Data(:, Index_Date), Data(:, Index_TTM));   
    clear Data_RF
    
    % Dividend Yield
    FileName = ['IndexDivYield19962021.txt'];
    Data_DY = load(FileName);
    clear FileName

    Data(:, Index_DY) = Data_DY(Data_DY(:, Index_Date)==Target_Date, end); 
    clear Data_DY
    
    % Implied Stock Price  
    K = Data(:, Index_K);
    OP = 0.5 * (Data(:, Index_OP_Bid) + Data(:, Index_OP_Ask));
    CPFlag = Data(:, Index_CPFlag);   
    S0 = Data(1, Index_S);
    TTM = Data(1, Index_TTM) / 365;                                                                           
    RF = Data(1, Index_RF);                                                    
    DY = Data(1, Index_DY);                                                   
     
    [AllK, ~, Index_AllK] = unique(K);                           
    Pair_OP_AllK = accumarray([Index_AllK CPFlag], ...
                              OP, ...
                              [length(AllK) 2], ...
                              @mean, ...
                              nan);                           
    Index_Drop = find(sum(isnan(Pair_OP_AllK), 2) > 0);
    if length(Index_Drop) > 0
        AllK(Index_Drop) = [];                                             
        Pair_OP_AllK(Index_Drop, :) = [];                                  
    else
    end
    clear Index_Drop   
    
    [~, Index_ATM] = min(abs(AllK / S0 - 1));
    K_ATM = AllK(Index_ATM);
    OP_C_ATM = Pair_OP_AllK(Index_ATM, 1);
    OP_P_ATM = Pair_OP_AllK(Index_ATM, 2);
    Data(:, Index_IS) = (OP_C_ATM - OP_P_ATM + K_ATM * exp(- RF * TTM)) / exp(- DY * TTM);
    clear Index_AllK AllK Pair_OP_AllK Index_ATM K_ATM OP_C_ATM OP_P_ATM

    Data(:, Index_IS_ADJ) = exp(- DY * TTM) * Data(:, Index_IS);  
    Data(:, Index_S_ADJ) = exp(- DY * TTM) * Data(:, Index_S); 
    clear K OP CPFlag S0 TTM RF DY

    % Data Filtering 
    Index = find((Data(:, Index_OP_Ask) > Data(:, Index_OP_Bid)) & ...
                 (Data(:, Index_OP_Bid) > 0)); 
    Data = Data(Index, :);                                                 
    clear Index

    Index = find(Data(:, Index_V) > 0);
    Data = Data(Index, :);                                                 
    clear Index      

    Index = find(Data(:, Index_OI) > 0);
    Data = Data(Index, :);                                                 
    clear Index      

    K = Data(:, Index_K);                                                                                             
    S0 = Data(:, Index_IS);      
    Index_C = find((Data(:, Index_CPFlag)==1) & ...
                   ((K ./ S0) >= 1) & ...
                   ((K ./ S0) <= 1.05));        
    Index_P = find((Data(:, Index_CPFlag)==2) & ...
                   ((K ./ S0) >= 0.9) & ...
                   ((K ./ S0) < 1));               
    Index = union(Index_C, Index_P);
    Data = Data(Index, :);                                                 
    clear Index Index_C Index_P K S0      
    
    OP = 0.5 * (Data(:, Index_OP_Bid) + Data(:, Index_OP_Ask));
    Index = find(OP >= (3 / 8));
    Data = Data(Index, :);                                                 
    clear Index OP      

    K = Data(:, Index_K);
    S0_ADJ = Data(:, Index_IS_ADJ);
    OP = 0.5 * (Data(:, Index_OP_Bid) + Data(:, Index_OP_Ask));
    TTM = Data(:, Index_TTM) / 365;                                        
    RF = Data(:, Index_RF);                                                    
    Index_C = find((Data(:, Index_CPFlag)==1) & ...
                   (OP <= S0_ADJ) & ...
                   (OP >= max(S0_ADJ - exp(- RF .* TTM) .* K, 0)));
    Index_P = find((Data(:, Index_CPFlag)==2) & ...
                   (OP <= (exp(- RF .* TTM) .* K)) & ...
                   (OP >= (exp(- RF .* TTM) .* K - S0_ADJ)));
    Data = [sortrows(Data(Index_C, :), Index_K); ...
            sortrows(Data(Index_P, :), Index_K)];                          
    clear K S0_ADJ OP TTM RF
    clear Index_C Index_P               
     
    % Put-Call Parity 
    Index = find(Data(:, Index_CPFlag)==2);                                    
    K = Data(Index, Index_K);
    S0_ADJ = Data(Index, Index_IS_ADJ);   
    TTM = Data(Index, Index_TTM) / 365;                                                                      
    RF = Data(Index, Index_RF);                                                     
    Data(Index, Index_CPFlag) = 1;                                                            
    Data(Index, Index_OP_Bid) = Data(Index, Index_OP_Bid) + S0_ADJ - K .* exp(- RF .* TTM);   
    Data(Index, Index_OP_Ask) = Data(Index, Index_OP_Ask) + S0_ADJ - K .* exp(- RF .* TTM);      
    Data = sortrows(Data, Index_K);                                        
    clear Index K S0_ADJ TTM RF    

    % Vertical Spread Condition
    Index = find(Data(:, Index_CPFlag)==1);                                
    K = Data(Index, Index_K);
    S0 = Data(Index(1), Index_IS);       
    CPFlag = Data(Index(1), Index_CPFlag);
    Index_Exist = CheckANAC_V(0.5 * (Data(Index, Index_OP_Bid) + Data(Index, Index_OP_Ask)), ...
                              K, S0, CPFlag);
    Index = Index(Index_Exist);                                                                   
    Data = Data(Index, :);                                                  
    clear Index K S0 CPFlag  
    clear Index_Exist   

    % Option-Implied Volatility
    Data(:, Index_IV_Bid) = nan;                                                                                      
    Data(:, Index_IV_Ask) = nan;                                                                                   
    Type = cell(size(Data, 1), 1);                                         
    Type(Data(:, Index_CPFlag)==1) = {'Call'};   
    Type(Data(:, Index_CPFlag)==2) = {'Put'};
    K = Data(:, Index_K);    
    S0_ADJ = Data(:, Index_IS_ADJ);   
    TTM = Data(:, Index_TTM) / 365;                                        
    RF = Data(:, Index_RF);                                                    
    Data(:, Index_IV_Bid) = blsimpv(S0_ADJ, K, RF, TTM, Data(:, Index_OP_Bid), [], [], [], Type);                                        
    Data(:, Index_IV_Ask) = blsimpv(S0_ADJ, K, RF, TTM, Data(:, Index_OP_Ask), [], [], [], Type);     
    clear Type K S0_ADJ TTM RF      
    
    % Output
    FileName = ['OP' num2str(Target_Date) '_TTM' num2str(Target_TTM) '.mat'];
    save(FileName, ...
         'Target_Date', 'Target_TTM', 'Data', ...
         'Index_Date', 'Index_TTM', 'Index_CPFlag', 'Index_K', 'Index_S', 'Index_OP_Bid', 'Index_OP_Ask', ...
         'Index_RF', 'Index_DY', 'Index_IS', 'Index_IS_ADJ', 'Index_S_ADJ', 'Index_IV_Bid', 'Index_IV_Ask');
    clear FileName

    clear Target_Date Target_TTM Data
    clear Index_ID Index_Date Index_TTM Index_CPFlag Index_K Index_S Index_F Index_OP_Bid Index_OP_Ask Index_OI Index_V
    clear Index_RF Index_DY Index_IS Index_IS_ADJ Index_S_ADJ Index_IV_Bid Index_IV_Ask
end
clear d